#include "StdAfx.h"
#include "NMFClass.h"

//Constructor
CNMFClass::CNMFClass(double** m_S, int m_C, int m_R, int m_NumberF, long MaxIte)
{
	//Determine if the original data is nonnegative or not
	for(int i=0; i<m_R; i++)
	{
		for(int j=0; j<m_C; j++)
		{
			if(m_S[i][j]<0)
			{
				CString str("Negative values of ");
				str += bf.douTostr(m_S[i][j]);
				str += CString(" at locations\r\nX: ");
				str += bf.douTostr(i);
				str += CString("\t, Y: ");
				str += bf.douTostr(j);
				::AfxMessageBox(str, MB_OK, 0);
				m_IsCreated = 0;
				return;
			}
		}
	}

	//Assign parameter to the NMF
	m_Original = m_S;
	m_Column = m_C;
	m_Row = m_R;
	m_NumberFactor = m_NumberF;
	m_MaxIteration = MaxIte;
	m_IsCreated = 1;
	
	//Scale the matrix to 0~1
	CNMFClass::NMFScale();

	//Initialize the NMF
	CNMFClass::NMFInitialize();

	return;
}

//Initialize the NMF
void CNMFClass::NMFInitialize(void)
{
	if(!m_IsCreated)
		return;
	//Create the basic and component
	m_Basic = (double**)malloc(sizeof(double*)*m_Row);
	m_Component = (double**)malloc(sizeof(double*)*m_NumberFactor);

	for(int i=0; i<m_Row; i++)
		m_Basic[i] = (double*)calloc(m_NumberFactor, sizeof(double));
	for(int i=0; i<m_NumberFactor; i++)
		m_Component[i] = (double*)calloc(m_Column, sizeof(double));

	//Create the approximation
	m_Approximation = (double**)malloc(sizeof(double*)*m_Row);
	for(int i=0; i<m_Row; i++)
		m_Approximation[i] = (double*)calloc(m_Column, sizeof(double));

	//Generate random doubles for the two matrix
	srand((unsigned)time(NULL));

	//Initialize the values of basic and component
	for(int i=0; i<m_Row; i++)
	{
		for(int j=0; j<m_NumberFactor; j++)
		{
			m_Basic[i][j] = (double)(rand()+1)/(double)(rand()+1);
			if(m_Basic[i][j]>1)
				m_Basic[i][j] = (double)1/m_Basic[i][j];
			//Move the matrix from 0~1 to 1~2
			m_Basic[i][j]+=1.0;
		}
	}
	for(int i=0; i<m_NumberFactor; i++)
	{
		for(int j=0; j<m_Column; j++)
		{
			m_Component[i][j] = (double)(rand()+1)/(double)(rand()+1);
			if(m_Component[i][j]>1)
				m_Component[i][j] = (double)1/m_Component[i][j];
			//Move the matrix from 0~1 to 1~2
			m_Component[i][j]+=1.0;
		}
	}

	//Determine the first approximation
	for(int i=0; i<m_Row; i++)
	{
		for(int j=0; j<m_Column; j++)
		{
			double NewApproximation = 0;
			for(int k=0; k<m_NumberFactor; k++)
				NewApproximation += m_Basic[i][k]*m_Component[k][j];
			m_Approximation[i][j] = NewApproximation;
		}
	}

	//Pre-assign previous cost to a large number and calculate the current cost
	m_PreviousCost = 999999999;
	m_CurrentCost = CNMFClass::NMFCalculateCost();

	return;
}

//Destroy the NMF
void CNMFClass::NMFDestroy(void)
{
	if(!m_IsCreated)
		return;
	//Free the memory of basics and approximation
	for(int i=0; i<m_Row; i++)
	{
		free(m_Basic[i]);
		free(m_Approximation[i]);
	}

	//Free the memory of components
	for(int i=0; i<m_NumberFactor; i++)
		free(m_Component[i]);

	m_IsCreated = 0;

	return;
}

//Factorizing the matrix
void CNMFClass::NMFFactorizing(void)
{
	if(!m_IsCreated)
		return;

	int m_Iteration = 0;
/*	double bmax=0;
	double bmin=999999;
	double cmax=0;
	double cmin=999999;
*/	//If unstable, train the system until the new cost is not longer decrease
	while(m_PreviousCost>m_CurrentCost && m_Iteration<m_MaxIteration)
	{
		//Update the basic
		for(int i=0; i<m_Row; i++)
		{
			for(int j=0; j<m_NumberFactor; j++)
			{
				//Calculate the first one
				double First = 0;
				for(int k=0; k<m_Column; k++)
					First += m_Original[i][k]*m_Component[j][k];
				//Calculate the second
				double Second = 0;
				for(int k=0; k<m_Column; k++)
					Second += m_Approximation[i][k]*m_Component[j][k];
				if(Second==0)
				{
/*					CString str("Second equal to 0 at location of i: ");
					str += bf.douTostr(i);
					str += CString(" j: ");
					str += bf.douTostr(j);
					str += CString("in basic update!!");
					::AfxMessageBox(str, MB_OK, 0);
*/					m_Basic[i][j]=m_Basic[i][j]*First;
				}
				else
					m_Basic[i][j] = m_Basic[i][j]*First/Second;
/*				if(m_Basic[i][j]>bmax)
					bmax=m_Basic[i][j];
				if(m_Basic[i][j]<bmin)
					bmin=m_Basic[i][j];
*/			}
		}
		//Update the component
		for(int i=0; i<m_NumberFactor; i++)
		{
			for(int j=0; j<m_Column; j++)
			{
				//Calculate the first one
				double First = 0;
				for(int k=0; k<m_NumberFactor; k++)
					First += m_Basic[i][k]*m_Original[k][j];
				//Calculate the second one
				double Second = 0;
				for(int k=0; k<m_NumberFactor; k++)
				{
					double m_Sum = 0;
					for(int m=0; m<m_Row; m++)
					{
						m_Sum += m_Basic[m][i]*m_Basic[m][k];
					}
					Second += m_Sum*m_Component[k][j];
				}
				if(Second==0)
				{
/*					CString str("Second equal to 0 at location of i: ");
					str += bf.douTostr(i);
					str += CString(" j: ");
					str += bf.douTostr(j);
					str += CString("in component update!!");
					::AfxMessageBox(str, MB_OK, 0);
*/					m_Component[i][j] = m_Component[i][j]*First;
				}
				else
					m_Component[i][j] = m_Component[i][j]*First/Second;
/*				if(m_Component[i][j]>cmax)
					cmax=m_Component[i][j];
				if(m_Component[i][j]<cmin)
					cmin=m_Component[i][j];
*/			}
		}
		//Update the approximation
//		double amin = 9999999;
//		double amax = 0;
		for(int i=0; i<m_Row; i++)
		{
			for(int j=0; j<m_Column; j++)
			{
				double NewApproximation = 0;
				for(int k=0; k<m_NumberFactor; k++)
					NewApproximation += m_Basic[i][k]*m_Component[k][j];
				m_Approximation[i][j] = NewApproximation;
/*				if(m_Approximation[i][j]<amin)
					amin=m_Approximation[i][j];
				if(m_Approximation[i][j]>amax)
					amax=m_Approximation[i][j];
*/			}
		}

		//Calculate the new cost
		m_PreviousCost = m_CurrentCost;
		m_CurrentCost = CNMFClass::NMFCalculateCost();
		::AfxMessageBox(bf.douTostr(m_CurrentCost), MB_OK, 0);
		m_Iteration ++;

/*		//Show the max and the min of the approximation
		::AfxMessageBox(CString("Basic max is: ")+bf.douTostr(bmax), MB_OK, 0);
		::AfxMessageBox(CString("Basic min is: ")+bf.douTostr(bmin), MB_OK, 0);
		::AfxMessageBox(CString("Component max is: ")+bf.douTostr(cmax), MB_OK, 0);
		::AfxMessageBox(CString("Component min is: ")+bf.douTostr(cmin), MB_OK, 0);
		::AfxMessageBox(CString("Approximation max is: ")+bf.douTostr(amax), MB_OK, 0);
		::AfxMessageBox(CString("Approximation min is: ")+bf.douTostr(amin), MB_OK, 0);
*/	}
	::AfxMessageBox(CString("Train finished at iteration: ")+bf.douTostr(m_Iteration), MB_OK, 0);
	return;
}

//Get the basic
double** CNMFClass::NMFGetBasic(int* m_C, int* m_R)
{
	*m_C = m_NumberFactor;
	*m_R = m_Row;
	return m_Basic;
}

//Get the component
double** CNMFClass::NMFGetComponent(int* m_C, int* m_R)
{
	*m_C = m_Column;
	*m_R = m_NumberFactor;
	return m_Component;
}

//Calculate the cost function
double CNMFClass::NMFCalculateCost(void)
{
	double TotalCost = 0;
	for(int i=0; i<m_Row; i++)
		for(int j=0; j<m_Column; j++)
		{
//			TotalCost += m_Original[i][j]*log(max(0.0000001, m_Original[i][j]/max(0.00000001,
//				m_Approximation[i][j])))-m_Original[i][j]+m_Approximation[i][j];
			TotalCost += (m_Original[i][j]-m_Approximation[i][j])*(m_Original[i][j]-m_Approximation[i][j]);
		}
	return TotalCost;
}

//Scale the matrix
void CNMFClass::NMFScale(void)
{
	//Determine the min and max
	m_Min = 999999999;
	m_Max = 0;
	for(int i=0; i<m_Row; i++)
	{
		for(int j=0; j<m_Column; j++)
		{
			if(m_Original[i][j]>m_Max)
				m_Max = m_Original[i][j];
			if(m_Original[i][j]<m_Min)
				m_Min = m_Original[i][j];
		}
	}

	if(m_Max == 0)
	{
		::AfxMessageBox(CString("Are you stupid!\r\nMax of matrix is 0"), MB_OK, 0);
		m_IsCreated = 0;
		return;
	}

	if(m_Max == m_Min)
	{
		::AfxMessageBox(CString("Are you stupid!\r\nThe elements of the matrix is the same"), MB_OK, 0);
		m_IsCreated = 0;
		return;
	}

	//Rescale the matrix to 1~2
	for(int i=0; i<m_Row; i++)
		for(int j=0; j<m_Column; j++)
			m_Original[i][j] = (m_Original[i][j]-m_Min)/(m_Max-m_Min)+1.0;
	
	return;
}
//Deconstructor
CNMFClass::~CNMFClass(void)
{
}
